setwd("...")
options("modelsummary_format_numeric_latex" = "plain")
library(dplyr)
library(tidyr)
library(modelsummary)
library(kableExtra)
library(ggplot2)
library(paletteer)
source("extra/func.R")
theme_set(my_theme())

data = read.csv("data/dataset_weekly.csv")

## Remove municipalities where there's only one observation?
data = subset(data, !(weeks_to_nc == 0 & week_since_str == 1))

### ====================================================================

## MODEL SPECS

# Formulae: effect of time (weeks since start) depending on presence of local orgs
f1 = "week_since_str + patr33_bin + CNT_bin + UGT_bin"
f2 = "week_since_str * patr33_bin + CNT_bin + UGT_bin"
f3 = "week_since_str * UGT_bin + patr33_bin + CNT_bin"
f4 = "week_since_str * CNT_bin + patr33_bin + UGT_bin"

# List
mtime = c(f1, f2, f3, f4)

# Controls
controls = c("pop1930_l", "izq1936", "analfab30_total", "clerics_all_l")
controls_f = paste(c(controls, "factor(prov_name)"), collapse = " + ")

### ====================================================================

## MAPS & CIA for TABLES

coef_recode = c(
  "patr33_bin" = "Employers' association",
  'CNT_bin' = 'CNT union',
  'UGT_bin' = 'UGT union',
  'pop1930_l' = 'Population (log)',
  'izq1936' = 'Leftist support 1936',
  'analfab30_total' = 'Illiteracy 1930',
  'clerics_all_l' = 'Number of clerics (log)',
  'week_since_str' = 'Weeks since war start',
  'week_since_str:patr33_bin' = "Weeks $\\times$ Employers' assoc",
  'week_since_str:CNT_bin' = "Weeks $\\times$ UGT union",
  'week_since_str:UGT_bin' = "Weeks $\\times$ CNT union"
)

gs = list(
  list("raw" = "nobs", "clean" = "$n$", "fmt" = 0),
  list("raw" = "r.squared", "clean" = "$R^2$", "fmt" = 2),
  list("raw" = "adj.r.squared", "clean" = "Adj. $R^2$", "fmt" = 2))

n = "+ p $<$ 0.1, * p $<$ 0.05, ** p $<$ 0.01, *** p $<$ 0.001. Municipality-week data, from July 17, 1936 to December 30, 1936. Control variables not shown: population 1930 (log), leftist share 1936, days of Republican control, illiteracy 1930, number of regular clerics (log), and number of secular priests (log). Province FE not shown. Only municipalities that were at least two weeks under Republican control."


### ====================================================================
## MODELS, main

lmw_all = lapply(mtime, function(x)
  lm(as.formula(paste0("weekly_ac_l ~", x, "+", controls_f)),
  data = data))

lmw_sec = lapply(mtime, function(x)
  lm(as.formula(paste0("weekly_ac_sec_l ~", x, "+", controls_f)),
  data = data))

lmw_reg = lapply(mtime, function(x)
  lm(as.formula(paste0("weekly_ac_reg_l ~", x, "+", controls_f)),
  data = data))

names(lmw_all) = paste0("(", 1:length(lmw_all), ")")
names(lmw_sec) = paste0("(", 1:length(lmw_sec), ")")
names(lmw_reg) = paste0("(", 1:length(lmw_reg), ")")

## Tables

# All clerics
modelsummary(
  models = lmw_all,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_rename = coef_recode,
  coef_omit = paste(c("prov_name", controls), collapse="|"),
  gof_map = gs,
  title = "Linear models on weekly anticlerical violence (logged), all clergy\\label{tab:lmw_all}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(lmw_all))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_weekly_all.tex")

# Secular
modelsummary(
  models = lmw_sec,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_rename = coef_recode,
  coef_omit = paste(c("prov_name", controls), collapse="|"),
  gof_map = gs,
  title = "Linear models on weekly anticlerical violence (logged), secular clergy\\label{tab:lmw_sec}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(lmw_sec))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_weekly_sec.tex")

# Regular
modelsummary(
  models = lmw_reg,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_rename = coef_recode,
  coef_omit = paste(c("prov_name", controls), collapse="|"),
  gof_map = gs,
  title = "Linear models on weekly anticlerical violence (logged), regular clergy\\label{tab:lmw_reg}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(lmw_reg))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_weekly_reg.tex")


## Predicted probabilities by week
pp = bind_rows(
  weekly_pp_int(modelnum = 2),
  weekly_pp_int(modelnum = 3),
  weekly_pp_int(modelnum = 4)) %>%
  mutate(org_lab = recode(org, "0" = "No local presence", "1" = "Local presence")) %>%
  mutate(model_lab = recode(model,
    "sind33_bin" = "General trade union",
    "patr33_bin" = "Employers' association",
    "CNT_bin" = "CNT",
    "UGT_bin" = "UGT"))

## Plots

# Main, only employers' assocs
p = ggplot(subset(pp, clergy == "all" & model == "patr33_bin"),
    aes(x = week_since_str, y = y, group = org_lab, fill = org_lab)) +
  geom_line(aes(color = org_lab), show.legend = FALSE) +
  geom_ribbon(aes(ymin = lwr90, ymax = upr90), alpha = 0.2) +
  geom_ribbon(aes(ymin = lwr95, ymax = upr95), alpha = 0.2) +
  labs(x = "Weeks since start of war", y = "Predicted clergy killings by week (log)\n") +
  theme(legend.title = element_blank(), legend.position = "bottom") +#c(0.9,0.9)) +
  scale_x_continuous(breaks = c(1, 5, 10, 15)) +
  scale_fill_manual(labels = c("Local employers' associations", "No local presence"),
    values = c("#FF0000", "#00A08A")) +
  scale_color_manual(values = c("#FF0000", "#00A08A"))
ggsave("output/predprob_weekly_patr33.pdf", height = 4, width = 5)

# Main model (all violence)
p = ggplot(subset(pp, clergy == "all"),
    aes(x = week_since_str, y = y, group = org_lab, color = org_lab)) +
  geom_line(position = position_dodge(1/5)) +
  geom_point(position = position_dodge(1/5)) +
  geom_errorbar(aes(ymin = lwr90, ymax = upr90), width = 0, linewidth = 1.1,
    position = position_dodge(1/5)) +
  geom_errorbar(aes(ymin = lwr95, ymax = upr95), width = 0, position = position_dodge(1/5)) +
  facet_wrap(~ model_lab, ncol = 4) +
  labs(x = "Weeks since start of war", y = "Predicted clergy killings by week (log)\n") +
  theme(legend.title = element_blank(), legend.position = "bottom") +#c(0.9,0.9)) +
  scale_x_continuous(breaks = c(1, 5, 10, 15)) +
  scale_color_paletteer_d(`"wesanderson::Darjeeling1"`)
ggsave("output/predprob_weekly.pdf", height = 4, width = 8)